home *** CD-ROM | disk | FTP | other *** search
/ The CICA Windows Explosion! / The CICA Windows Explosion! - Disc 2.iso / programr / dpmigcc5.zip / RSX / SOURCE / FPU-EMU / POLY_L2.C < prev    next >
C/C++ Source or Header  |  1994-05-27  |  8KB  |  290 lines
  1. /*---------------------------------------------------------------------------+
  2.  |  poly_l2.c                                                                |
  3.  |                                                                           |
  4.  | Compute the base 2 log of a FPU_REG, using a polynomial approximation.    |
  5.  |                                                                           |
  6.  | Copyright (C) 1992,1993                                                   |
  7.  |                       W. Metzenthen, 22 Parker St, Ormond, Vic 3163,      |
  8.  |                       Australia.  E-mail   billm@vaxc.cc.monash.edu.au    |
  9.  |                                                                           |
  10.  |                                                                           |
  11.  +---------------------------------------------------------------------------*/
  12.  
  13.  
  14. #include "exception.h"
  15. #include "reg_constant.h"
  16. #include "fpu_emu.h"
  17. #include "control_w.h"
  18.  
  19.  
  20.  
  21. #define    HIPOWER    9
  22. static unsigned short const    lterms[HIPOWER][4] =
  23.     {
  24.     /* Ideal computation with these coeffs gives about
  25.        64.6 bit rel accuracy. */
  26.     { 0xe177, 0xb82f, 0x7652, 0x7154 },
  27.     { 0xee0f, 0xe80f, 0x2770, 0x7b1c },
  28.     { 0x0fc0, 0xbe87, 0xb143, 0x49dd },
  29.     { 0x78b9, 0xdadd, 0xec54, 0x34c2 },
  30.     { 0x003a, 0x5de9, 0x628b, 0x2909 },
  31.     { 0x5588, 0xed16, 0x4abf, 0x2193 },
  32.     { 0xb461, 0x85f7, 0x347a, 0x1c6a },
  33.     { 0x0975, 0x87b3, 0xd5bf, 0x1876 },
  34.     { 0xe85c, 0xcec9, 0x84e7, 0x187d }
  35.     };
  36.  
  37.  
  38.  
  39.  
  40. /*--- poly_l2() -------------------------------------------------------------+
  41.  |   Base 2 logarithm by a polynomial approximation.                         |
  42.  +---------------------------------------------------------------------------*/
  43. void    poly_l2(FPU_REG const *arg, FPU_REG *result)
  44. {
  45.   short          exponent;
  46.   char          zero;        /* flag for an Xx == 0 */
  47.   unsigned short  bits, shift;
  48.   unsigned long long       Xsq;
  49.   FPU_REG      accum, denom, num, Xx;
  50.  
  51.  
  52.   exponent = arg->exp - EXP_BIAS;
  53.  
  54.   accum.tag = TW_Valid;    /* set the tags to Valid */
  55.  
  56.   if ( arg->sigh > (unsigned)0xb504f334 )
  57.     {
  58.       /* This is good enough for the computation of the polynomial
  59.      sum, but actually results in a loss of precision for
  60.      the computation of Xx. This will matter only if exponent
  61.      becomes zero. */
  62.       exponent++;
  63.       accum.sign = 1;    /* sign to negative */
  64.       num.exp = EXP_BIAS;  /* needed to prevent errors in div routine */
  65.       reg_u_div(&CONST_1, arg, &num, FULL_PRECISION);
  66.     }
  67.   else
  68.     {
  69.       accum.sign = 0;    /* set the sign to positive */
  70.       num.sigl = arg->sigl;        /* copy the mantissa */
  71.       num.sigh = arg->sigh;
  72.     }
  73.  
  74.  
  75.   /* shift num left, lose the ms bit */
  76.   num.sigh <<= 1;
  77.   if ( num.sigl & 0x80000000 ) num.sigh |= 1;
  78.   num.sigl <<= 1;
  79.  
  80.   denom.sigl = num.sigl;
  81.   denom.sigh = num.sigh;
  82.   poly_div4(&significand(&denom));
  83.   denom.sigh += 0x80000000;            /* set the msb */
  84.   Xx.exp = EXP_BIAS;  /* needed to prevent errors in div routine */
  85.   reg_u_div(&num, &denom, &Xx, FULL_PRECISION);
  86.  
  87.   zero = !(Xx.sigh | Xx.sigl);
  88.   
  89.   mul64(&significand(&Xx), &significand(&Xx), &Xsq);
  90.   poly_div16(&Xsq);
  91.  
  92.   accum.exp = -1;        /* exponent of accum */
  93.  
  94.   /* Do the basic fixed point polynomial evaluation */
  95.   polynomial((unsigned *)&accum.sigl, (unsigned *)&Xsq, lterms, HIPOWER-1);
  96.  
  97.   if ( !exponent )
  98.     {
  99.       /* If the exponent is zero, then we would lose precision by
  100.      sticking to fixed point computation here */
  101.       /* We need to re-compute Xx because of loss of precision. */
  102.       FPU_REG   lXx;
  103.       char    sign;
  104.       
  105.       sign = accum.sign;
  106.       accum.sign = 0;
  107.  
  108.       /* make accum compatible and normalize */
  109.       accum.exp = EXP_BIAS + accum.exp;
  110.       normalize(&accum);
  111.  
  112.       if ( zero )
  113.     {
  114.       reg_move(&CONST_Z, result);
  115.     }
  116.       else
  117.     {
  118.       /* we need to re-compute lXx to better accuracy */
  119.       num.tag = TW_Valid;        /* set the tags to Vaild */
  120.       num.sign = 0;        /* set the sign to positive */
  121.       num.exp = EXP_BIAS - 1;
  122.       if ( sign )
  123.         {
  124.           /* The argument is of the form 1-x */
  125.           /* Use  1-1/(1-x) = x/(1-x) */
  126.           significand(&num) = - significand(arg);
  127.           normalize(&num);
  128.           reg_div(&num, arg, &num, FULL_PRECISION);
  129.         }
  130.       else
  131.         {
  132.           normalize(&num);
  133.         }
  134.  
  135.       denom.tag = TW_Valid;    /* set the tags to Valid */
  136.       denom.sign = SIGN_POS;    /* set the sign to positive */
  137.       denom.exp = EXP_BIAS;
  138.       
  139.       reg_div(&num, &denom, &lXx, FULL_PRECISION);
  140.  
  141.       reg_u_mul(&lXx, &accum, &accum, FULL_PRECISION);
  142.  
  143.       reg_u_add(&lXx, &accum, result, FULL_PRECISION);
  144.       
  145.       normalize(result);
  146.     }
  147.  
  148.       result->sign = sign;
  149.       return;
  150.     }
  151.  
  152.   mul64(&significand(&accum),
  153.     &significand(&Xx), &significand(&accum));
  154.  
  155.   significand(&accum) += significand(&Xx);
  156.  
  157.   if ( Xx.sigh > accum.sigh )
  158.     {
  159.       /* There was an overflow */
  160.  
  161.       poly_div2(&significand(&accum));
  162.       accum.sigh |= 0x80000000;
  163.       accum.exp++;
  164.     }
  165.  
  166.   /* When we add the exponent to the accum result later, we will
  167.      require that their signs are the same. Here we ensure that
  168.      this is so. */
  169.   if ( exponent && ((exponent < 0) ^ (accum.sign)) )
  170.     {
  171.       /* signs are different */
  172.  
  173.       accum.sign = !accum.sign;
  174.  
  175.       /* An exceptional case is when accum is zero */
  176.       if ( accum.sigl | accum.sigh )
  177.     {
  178.       /* find 1-accum */
  179.       /* Shift to get exponent == 0 */
  180.       if ( accum.exp < 0 )
  181.         {
  182.           poly_div2(&significand(&accum));
  183.           accum.exp++;
  184.         }
  185.       /* Just negate, but throw away the sign */
  186.       significand(&accum) = - significand(&accum);
  187.       if ( exponent < 0 )
  188.         exponent++;
  189.       else
  190.         exponent--;
  191.     }
  192.     }
  193.  
  194.   shift = exponent >= 0 ? exponent : -exponent ;
  195.   bits = 0;
  196.   if ( shift )
  197.     {
  198.       if ( accum.exp )
  199.     {
  200.       accum.exp++;
  201.       poly_div2(&significand(&accum));
  202.     }
  203.       while ( shift )
  204.     {
  205.       poly_div2(&significand(&accum));
  206.       if ( shift & 1)
  207.         accum.sigh |= 0x80000000;
  208.       shift >>= 1;
  209.       bits++;
  210.     }
  211.     }
  212.  
  213.   /* Convert to 64 bit signed-compatible */
  214.   accum.exp += bits + EXP_BIAS - 1;
  215.  
  216.   reg_move(&accum, result);
  217.   normalize(result);
  218.  
  219.   return;
  220. }
  221.  
  222.  
  223. /*--- poly_l2p1() -----------------------------------------------------------+
  224.  |   Base 2 logarithm by a polynomial approximation.                         |
  225.  |   log2(x+1)                                                               |
  226.  +---------------------------------------------------------------------------*/
  227. int    poly_l2p1(FPU_REG const *arg, FPU_REG *result)
  228. {
  229.   char        sign = 0;
  230.   unsigned long long     Xsq;
  231.   FPU_REG          arg_pl1, denom, accum, local_arg, poly_arg;
  232.  
  233.  
  234.   sign = arg->sign;
  235.  
  236.   reg_add(arg, &CONST_1, &arg_pl1, FULL_PRECISION);
  237.  
  238.   if ( (arg_pl1.sign) | (arg_pl1.tag) )
  239.     {            /* We need a valid positive number! */
  240.       return 1;
  241.     }
  242.  
  243.   reg_add(&CONST_1, &arg_pl1, &denom, FULL_PRECISION);
  244.   reg_div(arg, &denom, &local_arg, FULL_PRECISION);
  245.   local_arg.sign = 0;    /* Make the sign positive */
  246.  
  247.   /* Now we need to check that  |local_arg| is less than
  248.      3-2*sqrt(2) = 0.17157.. = .0xafb0ccc0 * 2^-2 */
  249.  
  250.   if ( local_arg.exp >= EXP_BIAS - 3 )
  251.     {
  252.       if ( (local_arg.exp > EXP_BIAS - 3) ||
  253.       (local_arg.sigh > (unsigned)0xafb0ccc0) )
  254.     {
  255.       /* The argument is large */
  256.       poly_l2(&arg_pl1, result); return 0;
  257.     }
  258.     }
  259.  
  260.   /* Make a copy of local_arg */
  261.   reg_move(&local_arg, &poly_arg);
  262.  
  263.   /* Get poly_arg bits aligned as required */
  264.   shrx((unsigned *)&(poly_arg.sigl), -(poly_arg.exp - EXP_BIAS + 3));
  265.  
  266.   mul64(&significand(&poly_arg), &significand(&poly_arg), &Xsq);
  267.   poly_div16(&Xsq);
  268.  
  269.   /* Do the basic fixed point polynomial evaluation */
  270.   polynomial(&(accum.sigl), (unsigned *)&Xsq, lterms, HIPOWER-1);
  271.  
  272.   accum.tag = TW_Valid;    /* set the tags to Valid */
  273.   accum.sign = SIGN_POS;    /* and make accum positive */
  274.  
  275.   /* make accum compatible and normalize */
  276.   accum.exp = EXP_BIAS - 1;
  277.   normalize(&accum);
  278.  
  279.   reg_u_mul(&local_arg, &accum, &accum, FULL_PRECISION);
  280.  
  281.   reg_u_add(&local_arg, &accum, result, FULL_PRECISION);
  282.  
  283.   /* Multiply the result by 2 */
  284.   result->exp++;
  285.  
  286.   result->sign = sign;
  287.   
  288.   return 0;
  289. }
  290.